Gene Expression Array Analyses Predict Proto-Oncogene Expression During Perineural Invasion in Pancreatic Ductal Adenocarcinoma

Background/Aims: Pancreatic ductal adenocarcinoma is the tumor type with the highest incidence of perineural invasion. This study tries to identify the differentially expressed genes regulated between pancreatic ductal adenocarcinoma tissues with perineural invasion and without perineural invasion. Materials and Methods: The GSE102238 profile was downloaded. Gene function and pathway analysis were subsequently conducted. A protein–protein interaction network was constructed to search for hub genes. Both univariate Cox analysis and multivariate Cox analysis were calculated to identify prognostic factors. Quantitative real-time polymerase chain reaction (RT-PCR) and overall survival analysis of hub genes were used to verify. Results: Our study identified 242 differentially expressed genes including 68 upregulated differentially expressed genes and 174 downregulated differentially expressed genes, which were involved in important functions and pathways. Nine relevant core genes using protein–protein interaction analysis as well as nestin (NES)/vascular endothlial growth factor (VEGF) signaling pathway which is highly related to the pathological process of perineural invasion in pancreatic ductal adenocarcinoma were also discovered. The differentiation was identified as an independent prognostic factor (P < .05) after multivariate Cox analysis. Three upregulated genes (JUP, CALM1, and NES) and 6 downregulated genes (EPHA2, ARF1, ORM2, TERT, IL18, and CXCL3) were validated by quantitative RT-PCR and they all had markedly worse overall survival (P < .05). Conclusion: This analysis showed that 9 core genes including JUP, CALM1, NES, EPHA2, ARF1, ORM2, TERT, IL18, and CXCL3, as well as NES/VEGF signaling pathway, have a relationship with the development process of perineural invasion in pancreatic ductal adenocarcinoma. Cox analysis and overall survival analysis suggested differentiation as an independent prognostic factor and key roles for these 9 hub genes in perineural invasion prognosis in pancreatic ductal adenocarcinoma.


INTRODUCTION
Pancreatic cancer (PC) has become one of the most threatening tumors to human beings because of its difficult early diagnosis, strong invasion, low resection rate, high mortality, and short survival time.2][3] Even though the basic and clinical research on PC has made progress in recent years, the prognosis of PC is still grim, as well as the overall survival (OS) rate of 5 years is less than 5%. 4 Studies have shown that PC has a high incidence of perineural invasion (PNI), even up to 100%.This may be a way of metastasis of PC and an important cause of the recurrence of PC.Perineural invasion-positive patients are always associated with poor prognosis and low survival rate. 5,6rineural invasion means that cancer cells invade the adventitia and perineurium and even reach the neurointima and Schwann cells and neurons closely associated with it. 7Pancreatic ductal adenocarcinoma (PDAC) is the tumor type with the highest incidence of PNI, but the tumor size is not necessarily related to the occurrence of PNI, even if cancer can be seen only under the microscope, PNI can still occur.The occurrence of PNI is related to the neurophilicity of PC cells and the close anatomical location of the pancreas and nerve plexus.The distribution of the nerve plexus makes good contact between cancer cells and nerves.[10][11][12] Perineural invasion has been considered as an extension of lymphatic metastasis for many years because of the presence of lymphatic vessels in the adventitia.In recent years, it was found that the lymphatic vessels did not penetrate the epineurium, so the relationship between PNI and lymphatic metastasis was excluded. 13,14e molecular mechanism of PNI in PC is extremely complex, and most studies have focused on the interaction between nerve and tumor cells, and few have paid attention to the changes in the tumor matrix and microenvironment.The invasion of PC cells into peripheral nerves not only provides a pathway for the metastasis of PC but also leads to neural remodeling and changes in the neural environment, thus profoundly affecting the microenvironment of PC. 15,16 Through the expression profile GSE102238 from the Gene Expression Omnibus database (GEO) and bioinformatics analysis, this study was undertaken to elucidate the differentially expressed genes (DEGs) between 50 pairs of PDAC tissues and matched non-tumor tissues, of which 28 pairs were diagnosed as PNI by experienced pathologists, to explore the downstream molecules associated with PNI gene.By doing this, we hope to provide a novel understanding of the molecular basis of the effect of PNI on the microenvironment of PC.

MATERIALS AND METHODS
The study was approved by the ethics committee of the Chinese PLA General Hospital (approval no.2021PS213).

Gene Expression Microarray Data
The Gene expression profile GSE102238 was downloaded from the GEO (GEO, www.n cbi.n lm.ni h.gov /geo/ ).GSE102238 was based on Agilent-052909 CBC_lncRNAm-RNA_V3 platform.The GSE102238 dataset contained PDAC tissue from patients diagnosed with PNI (n = 28) and PDAC tissue from patients without PNI (n = 22).

Differentially Expressed Genes in Pancreatic Ductal
Adenocarcinoma Tissue Raw data were TXT files, which were assessed with GEO2R, which compares original submitter-processed data tables utilizing the GEO query and limma R packages (Bioconductor software).The volcano maps and box plots were completed using ggplot2 package R software and GEO2R online tools to illustrate the differential appearance.The box plots of the selected sample gene expression data before and after normalization are shown in Figure 1.Log Fold Change (FC) >0.5 or log FC <−0.5 as well as P < .05were cutoffs to determine DEGs between PDAC tissue samples from patients diagnosed with PNI and PDAC tissue specimens from patients without PNI.

Gene Ontology and Kyoto Encyclopedia of Genes and Genome Analyses
Cytoscape v3.4.0 (www.cytoscape.org) and ClueGO v2.33 were applied for identifying Gene Ontology (GO) and pathway analysis.Gene Ontology analysis including biological process (BP), molecular function (MF) and cellular component (CC), and Kyoto Encyclopedia of Genes and Genome (KEGG) analysis could determine the distinguishable genes' expression patterns and roles among these differentially expressed mRNAs.Overrepresented pathways showing P < .05were deemed statistically significant.

Gene Interaction Network Generation
Multiple DEGs in this work might be PNI associated and could be involved in PNI progression in the PDAC microenvironment.First, DEGs were entered in the Search Tool for the Retrieval of Interacting Genes (STRING) database (http: //www .string-db .org/), which yielded an interaction network (combined score above 0.4).Next, protein-protein interaction (PPI) networks comprising PNI-related genes of humans were obtained with Cytoscape v3.4.0.Core gene distribution in this network was generated with NetworkAnalyzer.Then, the Molecular Complex Detection (MCODE) plugin of Cytoscape was used for screening the network's modules.

Cox Regression Analysis
Univariate Cox analysis of OS with the "survival" package in R was conducted to identify prognostic factors The differentiation is an independent prognostic factor and key roles for these 9 prognostic genes in PNI prognosis in PDAC were validated.

Main Points
including clinical features with significant prognostic value, and P < .05 was considered to be statistically significant.The independent prognostic factors were identified by multivariate Cox regression analysis with the "survival" package in R. The clinical data were gathered from the GEO database and utilized to screen for the independent prognostic factor.The variables including age, gender, sizes, tumor node metastasis classification (TNM) stages, differentiation, localization of tumor, vessel invasion, and hub genes were included.P-values in univariate Cox analysis were calculated with log-rank method.If the previously mentioned variables with P < .05 in univariate Cox analysis, they would be included in multivariate Cox analysis, and the variables with P < .05 in multivariate Cox analysis could be regarded as statistically significance.Consequently, the differentiation was identified as an independent prognostic factor (P < .05).

Quantitative Real-Time Polymerase Chain Reaction
Total RNA was extracted by using Trizol Reagent (Ambion, Austin, TX, USA) from the 3 pairs of PDAC tissues with and without PNI of 6 patients with pathologically confirmed PDAC.The cDNA synthesis was performed using Reverse Transcription Kit (VAZYME, Nanjing, China).Quantitative RT-PCR was then carried out with the RT-PCR Kit (VAZYME).
Following the manufacturer's instructions, the thermal program underwent 10 minutes at 95°C, then 40 cycles of amplifications, 15 seconds at 95°C for denaturation, 60 seconds at 60°C for annealing, and 15 seconds at 95°C for the extension.The RT-PCR were carried out on QuantStudio 6 Thermal Cycler (ABI, Houston, TX, USA) using SYBR Green PCR Master Mix (VAZYME).An internal control including GAPDH, and the sequence of all the primers used have been described in Supplementary Table 1.Each sample was analyzed and calculated in triplicate.The 2 −ΔΔCt method was used to calculate the relative quantification of the hub genes.

Survival Analysis
In the OS assessment, 50 patients were assigned to the low and high groups based on the median expression levels of various hub genes.Then, Kaplan-Meier survival curves were generated with the "survival" package in R to compare the differences reported in the OS.

Differentially Expressed Genes in Pancreatic Ductal Adenocarcinoma Tissue
We downloaded the gene profile GSE102238 from the GEO and used GEO2R algorithm to confirm DEGs in PDAC tissue from patients diagnosed with PNI compared with PDAC tissue from patients without PNI, which were shown in the volcano plot (Figure 2).Using the cutoff criteria, 242 DEGs including 68 upregulated DEGs and 174 downregulated DEGs screened in PDAC tissue from patients diagnosed with PNI compared with PDAC tissue from patients without PNI were discovered based on the whole expression profile, top 10 DEGs were listed in Table 1 and a complete differential gene expression table was included as a supplementary file (Supplementary Table 2).

Gene Ontology Annotation Analysis of Differentially Expressed Genes
Functional analysis of the 242 DEGs was revealed using the Cytoscape software.Target genes were annotated to the GO pathway, which significantly enriched in the regulation of protein complex disassembly, regulation of blood vessel endothelial cell migration, positive regulation of smooth muscle cell proliferation, positive regulation of supramolecular fiber organization, positive regulation of actin filament polymerization and other biological processes (Figure 3).For MF, the DEGs were enriched in chemokine activity, methylated histone binding, transmembrane receptor protein tyrosine kinase activity, transmembrane-ephrin receptor activity, phosphoric diester hydrolase activity, phospholipase C activity, positive regulation of oxidoreductase activity, positive regulation of monooxygenase activity, and others.In addition,  GO CC analysis also showed that the DEGs were significantly enriched in the cortical actin cytoskeleton, specific granule, specific granule lumen, myelin sheath, calcium channel complex, and others.

Kyoto Encyclopedia of Genes and Genome Enrichment Analysis of Differentially Expressed Genes
Target genes were annotated to the KEGG pathway, which enriched in NF-κB signaling pathway, biosynthesis of the N-glycan precursor (dolichol lipid-linked oligosaccharide, LLO) and transfer to a nascent protein, uptake and actions of bacterial toxins, stimuli-sensing channels, RAB GEFs exchange GTP for GDP on RABs, Rab regulation of trafficking, signaling by VEGF, VEGFA-VEGFR2 pathway, VEGFR2mediated vascular permeability, signaling by ERBB2, PI3K/ AKT signaling in cancer, downregulation of ERBB2 signaling, opioid signaling, G-protein-mediated events, PLC betamediated events, EPH-ephrin signaling, EPHA-mediated growth cone collapse, EPH-ephrin-mediated repulsion of cells, and others.4).The 10 main high-degree hub nodes encompassed EPHA2, ABL1, NES, TERT, AGT, CALM1, ARF1, RASA1, EPHA4, and IL18.Of the abovementioned genes, EPHA2 had the highest node degree of 10.Table 3 shows core genes and their respective degrees.Figure 5 depicts the overall prospect of the complex regulatory relationship among the detected core genes.The data points and the respective points on the graph were highly correlated (coefficient approximating 0.981).An R 2 of 0.986 was obtained, indicating the linearity of the model.Next, MCODE was utilized for screening modules in the gene interaction network, 4 of which are depicted in Figure 6.The first module comprising CXCL3, HTR1B, and GNAI1 had a score of 3, with 3 nodes and 3 edges.The second module also had 3 nodes and 3 edges and encompassed FCAR, SLC44A2, and TNFRSF1B, with a score of 2. The third module (JUP, ORM2, and ERP44) has a score of 2, and included 3 nodes and 3 edges.Finally, the fourth module (SRSF10, SNRPA, HNRNPDL, TERT, BMI1, and NES) had a score of 2, with 6 nodes and 7 edges.

Validation of Quantitative Real-Time Polymerase Chain Reaction
In addition to validating the bioinformatic analysis results, quantitative RT-PCR was used to quantify parts of explored genes, including 3 upregulated genes (JUP, CALM1, and NES) and 6 downregulated genes (EPHA2, ARF1, ORM2, TERT, IL18, and CXCL3).As shown in Figure 7, the gene expression patterns of JUP, CALM1, NES, EPHA2, ARF1, ORM2, TERT, IL18, and CXCL3 detected by quantitative RT-PCR significantly were accorded with the corresponding gene alteration of microarray data (P < .05).

Survival Analysis
Cases in the GSE102238 dataset were assigned to 2 groups, based on the median expression levels of various hub genes and Kaplan-Meier survival analysis.For the hub genes JUP, CALM1, and NES, respectively, individuals with elevated gene expression had markedly worse OS (P < .05).For the hub genes EPHA2, ARF1, ORM2, TERT, IL18, and CXCL3, respectively, individuals with lower gene expression showed significantly worse OS (P < .05;matrix metalloproteinases, and other effectors; then, these genes undergo changes in the generation, and tumor cells subsequently invade the nerve tissue. 17,18hang et al's 19 study identified PNI-associated genes in PC cell lines and the "chemokine signaling pathway" was found to be associated with PNI, following KEGG pathway enrichment analysis and the construction of a PPI network from the identified DEGs.Furthermore, FGF2 was found to be associated with PNI.Li et al's 20 review showed gene alternations in human PDAC samples are also linked to PNI.For example, Ras homolog family member C was abundantly expressed in PNI tissues and was related to poor disease prognosis.Here, GSE102238's gene profiling data were obtained and analyzed by bioinformatics.As  shown earlier, 242 DEGs were identified in PDAC tissue specimens from patients diagnosed with PNI compared with PDAC tissue samples from patients without PNI.In addition, enrichment analysis and gene interaction network analyses were carried out for identifying biomarkers or key genes associated with cytogenetic pathways controlling PNI in PDAC.Univariate Cox and multivariate Cox regression analyses were used to identify the most significant prognostic factors and genes.Finally, the differentiation was identified as an independent prognostic factor (P < .05)and 3 upregulated genes (JUP, CALM1, and NES) and 6 downregulated genes (EPHA2

Figure 1 .
Figure 1.The box plots of the selected sample gene expression data before (A) and after (B) normalization in the GSE102238 dataset.Green bars: pancreatic ductal adenocarcinoma (PDAC) tissue samples from patients diagnosed with perineural invasion (PNI); purple bars: PDAC tissue samples from patients diagnosed without PNI.

Figure 2 .
Figure2.Volcano plot of all differentially expressed genes (DEGs) between pancreatic ductal adenocarcinoma (PDAC) tissue samples from patients diagnosed with perineural invasion (PNI) and without PNI in the GSE102238 dataset.Red dots: significantly upregulated genes; blue dots: significantly downregulated genes; black dots: nondifferentially expressed genes.P < .05 was considered statistically significant.

Figure 3 .
Figure 3. Gene Ontology (GO)-enrichment analysis of biological processes (A), molecular functions (B), and cellular components (C).The red star in GO terms means Term P < .05,and double red stars in GO terms mean Term P < .01.

Figure 4 .
Figure 4. Protein-protein interaction network of differentially expressed genes.

Figure 8 )Figure 6 .
Figure8).These findings suggested key roles for these genes in PNI prognosis in PDAC.

Figure 5 .
Figure 5.The distribution of core genes in the interaction network.The black node means the core gene.The red line means the fitted line and the blue line means the power law.The correlation between the data points and corresponding points on the line is approximately 0.981.The R-squared value is 0.916 giving a relatively high confidence that the underlying model is indeed linear.

Figure 7 .
Figure 7. Results of the quantitative real-time PCR.Perineural invasion (PNI) groups: mRNA levels of pancreatic ductal adenocarcinoma (PDAC) tissue samples from patients diagnosed with PNI; control groups: mRNA levels of PDAC tissue samples from patients diagnosed without PNI.Asterisk means P < .05.

Figure 8 .
Figure 8. Overall survival analyses 9 hub genes in data set GSE102238.

Table 1 .
The Top 10 Regulated Differentially Expressed Genes in Pancreatic Ductal Adenocarcinoma Tissues with Perineural Invasion and Without Perineural Invasion with P < .05 DEGs, differentially expressed genes; FC, fold change; PDAC, pancreatic ductal adenocarcinoma; PNI, perineural invasion.

Table 2 .
Core Pathways and Their Associated Genes Found GO, Gene Ontology.

Table 3 .
The Core Genes and Their Corresponding Degree

Table 1 .
It also promotes the specific growth tendency of axons and increases nerve-tumor contact.[23][24][25][26]Theabovementionedpathwaysmayhave critical functions in PNI development in PDAC., and drug sensitivity in PDAC.This could be explained by the fact that elevated amounts of tumor vessels increase the odds of cancer cells to enter the bloodstream.New tumor vessels or capillaries possess weak and leaky basement membranes, and cancer cells could penetrate them more readily compared with mature counterparts.Furthermore, cancer cells can directly invade the nerve and can also invade the nerve through the penetrating channels (such as blood vessels and reticular fibers).33Ourresultalso discovered that the VEGF signaling pathway was one of the most relevant pathways for PNI in PDAC.So nestin may mediate increased PNI in PDAC by raising tumor cells invading the nerve through the blood vessels and activating the VEGF signaling pathway.Nestin could act as a novel therapeutic target for PC via tumor angiogenesis.Author Contributions: Concept -S.Z., P.S.; Design -S.Z., P.S.; Supervision -P.S.; Resources -P.S.; Materials -S.Z., P.S.; Data Collection and/or Processing -S.Z., Z.X., J.W.; Analysis and/or Interpretation -S.Z., Z.X., J.W.; Literature Search -S.Z., Z.X.; Writing Manuscript -S.Z., Z.X.; Critical Review -P.S., J.W. Primers of Genes Validated by qRT-PCR 22Overexpressed CXCR4 in bile duct carcinoma and PC is closely related to PNI, lymphatic metastasis, TNM staging, and vessel invasion.CXCR4 promotes VEGF expression, the mitosis and proliferation of vascular endothelial cells, and the tumor to generate new vessels.SNRPA, HNRNPDL, TERT, BMI1, and NES were detected as core interaction genes, which might constitute targets for treating PNI in PDAC.Of these, TERT and NES corroborated STRING database findings.Nestin (NES) is a class VI intermediate filament protein, which is distributed in the cytoplasm and involved in the formation of the cytoskeleton.It was previously considered as a marker of neural stem cells.In recent years, it has been found that nestin is expressed in pancreatic stem cells and PC stem cells and is related to tumorigenesis, tumor angiogenesis, tumor metastasis, and prognosis of PC.27The neural invasion pathway of PC is usually through the direct destruction of the perineural membrane, the vascular invasion of the perineural membrane, and the destruction of the synaptic membrane of the nerve endings.28Inthisstudy,wefocused on analyzing the perineural invasiveness of PDAC based on nestin expression in cancer cells.Kawamoto and collaborators as well as other investigators revealed nestin amounts in tumor cells correlate with nerve invasion in PC.Ethics Committee Approval: The study was approved by the ethics committee of the Chinese PLA General Hospital (approval no.2021PS213, Beijing, China).Peer-review: Externally peer-reviewed.2011;11(10):695-707. [CrossRef] 4. Sivapalan L,Kocher HM,Ross-Adams H, Chelala C. Molecular profiling of ctDNA in pancreatic cancer: opportunities and challenges for clinical application.Pancreatology.2021;21(2):363-378.[CrossRef] 5. Schorn S, Demir IE, Haller B, et al.The influence of neural invasion on survival and tumor recurrence in pancreatic ductal adenocarcinoma -a systematic review and meta-analysis.Surg Oncol.2017;26(1):105-115.[CrossRef] 6. Alrawashdeh W, Jones R, Dumartin L, et al.Perineural invasion in pancreatic cancer: proteomic analysis and in vitro modelling.Mol Oncol.2019;13(5):1075-1091. [CrossRef] 7. Hirai I, Kimura W, Ozawa K, et al.Perineural invasion in pancreatic cancer.Pancreas.2002;24(1):15-25.[CrossRef] 8. Demir IE, Friess H, Ceyhan GO.Neural plasticity in pancreatitis and pancreatic cancer.Nat Rev Gastroenterol Hepatol.2015;12(11):649-659. [CrossRef] 9. Ceyhan GO, Demir IE, Altintas B, et al.Neural invasion in pancreatic cancer: a mutual tropism between neurons and cancer cells.between multiple clinicopathological features and nerve invasion in pancreatic cancer.Hepatobiliary Pancreat Dis Int.2013;12(5):546-551.[CrossRef] 15.Noelle J,Lei Z. Signaling in the microenvironment of pancreatic cancer: transmitting along the nerve.Pharmacol Ther.2019;200:126-134.[CrossRef] 23.Zaitseva L, Murray MY, Shafat MS, et al.Ibrutinib inhibits SDF1/ CXCR4 mediated migration in AML.Oncotarget.2014;5(20):9930-9938. [CrossRef] 24.Li J, Ma Q, Liu H, et al.Relationship between neural alteration and perineural invasion in pancreatic cancer patients with hyperglycemia.PLOS ONE.2011;6(2):e17385.[CrossRef] 25.Li X, Ma Q, Xu Q, et al.SDF-1/CXCR4 signaling induces pancreatic cancer cell invasion and epithelial-mesenchymal transition in vitro through non-canonical activation of Hedgehog pathway.31.Bao Y, Prescott J, Yuan C, et al.Leucocyte telomere length, genetic variants at the TERT gene region and risk of pancreatic cancer.Gut.2017;66(6):1116-1122. [CrossRef] 32.Faleiro I, Apolónio JD, Price AJ, et al.The tert hypermethylated oncologic region predicts recurrence and survival in pancreatic cancer.Future Oncol.2017;13(23):2045-2051.[CrossRef] 33.Campa D, Rizzato C, Stolzenberg-Solomon R, et al.Tert gene harbors multiple variants associated with pancreatic cancer susceptibility.Int J Cancer.2015;137(9):2175-2183.[CrossRef] 34.Yamahatsu K, Matsuda Y, Ishiwata T, Uchida E, Naito Z. Nestin as a novel therapeutic target for pancreatic cancer via tumor angiogenesis.Int J Oncol.2012;40(5):1345-1357. [CrossRef] Supplementary

Table 2 .
The Complete Regulated DEGs in PDAC Tissues with PNI and without PNI with P < .05(Continued)

Table 2 .
The Complete Regulated DEGs in PDAC Tissues with PNI and without PNI with P < .05(Continued)

Table 2 .
The Complete Regulated DEGs in PDAC Tissues with PNI and without PNI with P < .05(Continued)